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I.  INTRODUCTION 


Background 

Progress  in  obtaining  general  solutions  to  viscous  fluid  flow 
problems  has  been  very  slow  and  difficult.  It  was  over  230  years  ago 

that  Euler  (1755)  derived  the  equations  for  frictionless  fluid  flow.1 
It  was  another  70  years  before  Navier  (1827)  added  viscous  terms  to 
Euler’s  equations.  Soon  after  that  Stokes  (1845)  improved  Navier ’s  idea 
by  the  introduction  of  the  coefficient  of  viscosity,  n,  that  defined  the 
fluid  stress  as  proportional  to  the  rate-of -change  of  strain.  Thus,  the 
equations  governing  the  flow  of  a  viscous  fluid  (Navier-Stokes 
equations)  have  been  known  for  over  140  years.  However,  only  a  few 
special-case  solutions  to  these  equations  have  been  found  analytically. 
The  equations  are  simply  too  difficult  to  solve  analytically  for  general 
flow  problems.  Fortunately  the  recent  development  of  large,  efficient 
computers  has  opened  the  way  for  numerical  solutions  of  the  governing 
equations . 

This  section  will  show  how  the  steady,  two-dimensional  Navier- 
Stokes  equations  are  formulated  for  numerical  solutions  and  how  this 
solution  is  obtained.  Much  of  this  information  is  available  from  text 
books  and  other  sources  but  it  is  summarized  here  for  convenience  and 
completeness.  Also  included  is  a  discussion  of  the  necessity  for,  and 
the  application  of,  a  turbulence  model  for  estimating  the  effects  of  a 
turbulent  boundary  layer. 

The  remainder  of  the  report  describes  the  modifications  made  to  the 
turbulence  model  in  order  to  get  a  successful  solution,  and  a  discussion 
of  the  results  for  transonic  flow  over  a  thick  supercritical  airfoil. 

The  thick  supercritical  airfoil  was  chosen  for  this  study  because  it 
provided  the  most  severe  flow  conditions  the  author  could  find  for  which 
experimental  data  were  available. 
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Governing  Equations 


The  Navier-Stokes  equations  can  be  derived  by  applying  the 
principles  of  conservation  of  mass,  conservation  of  momentum,  and  con- 

2 

servation  of  energy  to  a  control  volume  .  Using  an  Euler ian  formulation 
and  assuming  an  arbitrary  control  volume  V  with  surface  area  A,  we  can 
express  the  conservation  of  mass  as: 

Rate  of  mass  loss  in  the  control  volume 
equals  net  flow  of  mass  out  through 
the  surface. 

-  h  UIp  av  =  //  P  5  .  dX  (i) 

Where , 

t  =  time 

p  =  density  (mass  per  unit  volume) 

5  =  total  velocity  vector. 


Similarly,  the  conservation  of  momentum  (or  Newton’s  second  law)  states: 
The  sum  of  the  external  forces  equals  the 
rate  of  change  of  momentum. 


Where, 


II  p  •  dX  =  ///  p 


DU 

Ct 


dV 


(2) 


P  =  a  dyadic  representing  the  normal  and  shear 
stresses  acting  on  the  control  volume.  All 
other  forms  of  external  forces  are  assumed 
to  be  negligible. 

=  the  substantial  derivative  which  represents 

the  sum  of  the  local  change  in  the  variable  plus  the 
convective  rate-of-change . 


DO  _  90  & 
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Finally,  the  conservation  of  energy  states: 

The  rate-of-change  of  energy  in  the  control 


2 


volume  equals  the  rate  of  heat  added  through 
the  surface  plus  the  rate  of  work  done  by 
external  forces. 

_  ■*  ■* 

III  dV  =  //  q-dA  +  fj  (P.U).dA  (3) 

Where, 

e  =  internal  energy  per  unit  mass  (specific 
internal  energy  plus  kinetic  energy,  all  other 
energy  terms  are  assumed  to  be  negligible) 

-» 

q  =  heat  flux  vector. 

This  formulation  assumes  there  is  no  energy  source  or  sink  within  the 
control  volume. 

Notice  that  equation  (2)  is  a  vector  equation  and  represents  the 
conservation  of  momentum  in  each  dimension  of  interest  (x  and  y  for  a 
two-dimensional  problem  in  Cartesian  coordinates) .  Also  at  the  begin¬ 
ning  of  this  section  it  was  stated  that  the  control  volume  was  of 
arbitrary  size;  however,  it  is  assumed  to  contain  enough  molecules  to 

3 

represent  a  continuum. 

It  is  possible  to  develop  a  numerical  solution  procedure  to  solve 

the  governing  equations  in  integral  form.  For  example,  Jameson*  uses 
Runge-Kutta  integration  techniques  over  finite-volumes  to  solve  the 
Euler  equations  (equations  1-3  without  the  viscous  stresses).  However, 
the  numerical  method  used  in  the  present  study  requires  a  differential 
form  of  the  equations.  This  is  accomplished  by  applying  Gauss’  diver¬ 
gence  theorem  to  the  integral  equation  and  making  an  argument  about  the 
arbitrary  size  of  the  control  volume  V.  Gauss’  divergence  theorem 
equates  surface  integrals  to  volume  integrals  and  is  written  as: 

//  $  .  dA  =  ///  (V  .  j)  dV.  (4) 

Where , 

$  =  any  tensor  quantity  that  is  first-order  continuous, 

V  =  the  divergence  operator. 

Applying  Gauss’  divergence  theorem  to  the  right-hand  side  of  equation 
(1)  yields; 
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///  (V./>U)  dv  ♦  ^  ;//  p  av  =  0.  (5) 

If  the  control  volume  is  not  a  function  of  time,  the  order  of  integra¬ 
tion  and  differentiation  may  be  reversed  and  equation  (5)  may  be  written 


///  [(V./>U)  dV  =  0.  (6) 

At  this  point  the  argument  is  made  that  if  the  control  volume  is  ar¬ 
bitrary,  the  integrand  itself  must  vanish  if  the  integral  is  to  be  zero. 
Therefore, 

^  +  V  .  pt}  =  0  (7) 


is  the  derivative  form  of  the  conservation  of  mass, 
procedure  equation  (2)  may  be  written  as; 

and  equation  (3)  may  be  written  as; 


P  Dt 


De  -  V.(?.U  *  q)  =  0 


Following  a  similar 


(8) 

(9) 


Equations  (8)  and  (9)  still  contain  substantial  derivative  terms. 
These  may  be  written  in  a  more  convenient  form  by  using  the  definition 
of  the  substantial  derivative  and  making  use  of  the  conservation  of  mass 
equation.  The  equations  then  become: 


Continuity 


^  +  V.pU  =  0 


at 


(10) 


Momentum 

+  V.(pUU  -  5)  =  0 


(11) 


Energy 


•+  -» 

V. (pVe  -  P.C  -  q)  =  0. 


(12) 
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Equations  (10) ,  (11)  and  (12)  are  the  classic  Navier-Stokes  equations 
for  compressible  fluid  flow  without  body  forces  or  energy  loss  or  gain 
within  the  domain  of  interest.  However,  the  system  of  equations  is  not 
closed  and  several  variables  need  clarification  before  practical  solu¬ 
tions  can  be  obtained. 

-► 

The  heat  flux  (q)  is  determined  from  Fourier’s  law  as; 
q  =  k  VT 

Where 

T  =  the  temperature 

fc  =  the  thermal  conductivity. 

The  value  of  fc  for  fluids  is  a  thermodynamic  property  that  varies  as  a 
function  of  temperature. 

The  most  complicated  variable  in  the  Navier-Stokes  equations  is  the 
■*  ■* 
stress  tensor  (?) .  In  a  two-dimensional  Cartesian  coordinate  system  P 
can  be  written  as 


r 

T 

XX 

xy 

T 

r 

yx 

yy 

(13) 


the  normal  stress  in  the  x  direction 

the  normal  stress  in  the  y  direction 

the  shear  stress  on  the  x  face  in  the  y  direction 

the  shear  stress  on  the  y  face  in  the  x  direction. 

It  is  at  this  point  that  Stokes’  work  becomes  very  important.  He  postu¬ 
lated  that  the  stress  is  a  linear  function  of  strain  rates,  that  a  fluid 
is  isotropic,  and  that  the  normal  stress  must  equal  the  hydrostatic 
pressure  when  the  fluid  is  at  rest.  Using  these  assumptions  and  defin¬ 
ing  two  additional  variables  fi  and  X,  assumed  to  be  properties  of  the 
fluid,  the  stresses  may  be  defined  as 


where , 
r 

xx 

r 

yy 

r 

xy 

r 

yx 
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where  p  is  the  hydrostatic  pressure,  ft  is  the  coefficient  of  shear  vis¬ 
cosity  and  is  a  thermodynamic  property  of  the  fluid,  usually  computed 

from  Sutherland’s  formula*.  X  is  much  more  controversial.  Fortunately 
the  value  for  X  has  little  effect  on  the  solutions  of  interest  in  this 
report.  The  most  commonly  accepted  method  of  determining  the  value  of  X 
for  air  is  to  again  follow  Stokes’  lead  and  set  X  =  -(2/3)  p.  This 

forces  the  trace  of  ?  to  be  the  hydrostatic  pressure,  p. 

The  velocity  vector  C  for  two-dimensional  Cartesian  coordinates  may 

be  defined  in  terms  of  u,  the  component  in  the  x  direction  and  v,  the 
component  in  the  y  direction.  Thus  equation  (11)  may  be  separated  into 
two  scalar  equations.  This  leaves  a  system  of  four  equations  with  six 
unknowns  (u,  v,  p,  T,  e,  p) .  Obviously  some  other  relationships  are 
required.  The  additional  relationships  may  be  determined  by  assuming 
the  fluid  of  interest  (air)  obeys  the  perfect  gas  law.  This  provides  a 
relationship  between  density,  pressure  and  temperature, 

p  =  pRT  (15) 

where  R  is  the  gas  constant,  a  known  quantity  for  a  given  gas.  The  per¬ 
fect  gas  assumption  also  provides  a  relationship  between  internal  energy 
and  temperature, 


where, 

e^  =  the  specific  internal  energy  per  unit  mass 
Cy  =  the  constant  volume  specific  heat. 


I 


i 


1 


1 
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While  in  reality  Cy  is  a  function  of  temperature,  it  is  essentially  con¬ 
stant  for  practical  transonic  applications.  The  total  internal  energy 
is  the  sum  of  the  specific  internal  energy  and  the  kinetic  energy. 

e  =  CyT  ♦  1/2  (u2  +  v2)  (16) 

This  provides  a  system  of  six  equations  and  six  unknowns.  Notice  that 
equations  (IS)  and  (16)  may  be  used  to  compute  pressure  and  temperature 
in  terms  of  the  other  flow  variables.  Thus  a  system  of  four  equations 
in  four  unknowns  (u,  v,  p,  e)  may  be  established.  These  four  equations 
are  repeated  here  as  they  would  be  written  for  a  two-dimensional 
Cartesian  coordinate  system. 

Continuity 

*1  +  0£v  o  (17a) 

0t  9x  0y  u 


Momentum 


0/7U  9  ,  2 

St  * 


t  )  +  -s-  (py u  -  r  )  =  0 
xx'  9yvr  yx' 


(17b) 


9  pv 
Ot 


*  Si*'”  -  V 


-  T  )  =  0 

yy 


(17c) 


Energy 


dpe  9  ,  .  0T. 

sr—  +  w—  (/>ue  -  ur  -  vr  -  ks— ) 
ot  ox  r  xx  yx  0x' 

-5~(pve  -  vr 

Oy  r  yy 


UT  -  kE)  =  0 

xy  Oy' 


(17d) 


It  is  convenient  to  define  the  primary  flow  variables  as  p,  pn,  py ,  and 
pe.  Temperature  and  pressure  are  inherent  in  the  equations,  but  can  be 
computed  from  the  four  primary  variables  using  equations  (16)  and  (16)  . 
It  is  also  convenient  to  express  equations  (17)  in  flux  vector 

form. 


9U  9F 
8t  +  97  + 


(18) 


I 

1 


I 
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Where  U  here  is  not  the  velocity  vector  but  is  the  vector  of  the  depend¬ 
ent  variables.  This  traditional  convention  is  retained  even  at  the  risk 
of  some  confusion.  F  and  G  are  flux  vectors. 


r  =  r  is  assumed  for  the  stress  terms . 
xy  yx 

Generalized  Transformation 

Flow  solutions  may  be  computed  by  numerically  solving  equation  (18) 
on  an  appropriate  Cartesian  computational  grid.  In  general,  however, 
the  body  about  which  the  flow  solution  is  desired  does  not  coincide  with 
the  grid  lines  of  a  Cartesian  grid.  As  a  result,  it  is  very  difficult 
to  obtain  adequate  grid  resolution  near  the  solid  body  surface,  and  the 
boundary  conditions  at  the  surface  are  difficult  to  enforce  accurately. 
Therefore,  a  body  conforming  computational  grid  is  usually  established 
and  the  necessary  transformations  from  the  Cartesian  to  the  computa¬ 
tional  coordinate  system  are  performed.  The  resulting  equation  is 


8 


similar  to  equation  (18),  but  tedious  mathematics  are  required  to  obtain 
a  generalized  transformation.  Use  of  matrix  representation  and  the 


notation  jjjj  =  will  be  used  to  help  simplify  the  following  develop¬ 


ment. 

The  grid  points  in  the  body  conforming  grid  are  usually  defined  in 
terms  of  locations  in  a  Cartesian  coordinate  system.  That  is,  if  f  and 
1 1  represent  the  transformed  coordinate  system,  the  computational  grid  is 
defined  as; 

(  =  £(x,y) 

V  =  >?(x,y) 


The  total  differentials 
d{  =  (x  dx  + 

dty  =  i?x  dx  + 

Or  in  matrix  form; 


for  {  and  tf  are  defined  by; 
f  dy 

Vy  dy 


V 

1 

X 

X 

v</» 

1 _ 

dx 

dq 

V  V 
_'x  'y_ 

dy 
-  — 

(19) 


The  reverse  transformation  is: 
x  =  x(£,I7) 

y  *  y (C»*7) 

and  the  matrix  looks  like: 


dx 

x-  X 

d£ 

dy 

l 

St- 

X 

w 

X 

_ 1 

d  V 

(20) 


Multiplying  (20)  by  the  inverse  of  the  partial  derivative  matrix  and 
comparing  the  results  with  (19)  provides: 


V 

dx 

’{  "i 

-1 

dx 

d*? 

\  *y 

dy 

r{  y-» 

—  -am 

dy 

I 


I 


I 


I 
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That  is,  the  matrix  of  the  partial  derivatives  of  the  body  conforming 
grid  is  equal  to  the  inverse  of  the  matrix  of  partial  derivatives  of  the 
Cartesian  coordinates.  Linear  algebra  techniques  provide: 


If  the  determinant,  called  the  Jacobian  (J)  of  the  transformation, 


is  non-zero,  a  one-to-one  mapping  exists  such  that: 


J  ut  *v  -  V’f  * 


(«{G  -  y{P), 


(-Py  *  +  Gx  #  -  Cx.  +  Fy,  )  =  0 

V  *Vi  V(  iv  ivJ 


(23) 


If  the  functions  x  =  x(£,r))  and  y  =  y(£,v)  are  continuous  to  at 


least  the  second  derivative,  x  *  =  x,.  and 

vi  iv 


y  ,  -  y,  and  the  last  terms 

vi  iv 


in  parenthesis  in  equation  (23)  vanish.  The  finite  difference  repre¬ 
sentation  (discussed  below)  of  these  metrics  must  also  vanish  or  errors 

will  be  introduced  into  the  finite  difference  solutions5. 

Redefine  the  variables  in  equation  (23)  such  that: 

U  =  J  U 

F  =  y  F  -  x  G 

V  V 

G  =  x^G  -  y^F 

and  the  transformed  equation  of  motion  is: 

=  0 


9U 

dt 


8F  8G 

+  -s-r-  + 


(24) 


H  a* 

This  is  the  strong  conservation  form  of  the  transformed  governing  equa¬ 
tions.  Conservation  form  is  defined  by  Ref.  5  as  a  form  where  the 
derivatives  of  the  dependent  variables  have  only  constant  coefficients. 
Strong  conservation  form  means  there  are  no  source  terms  (right  hand 
side  =  0) .  Conservation  is  very  important  for  correctly  computing 
strong  flow  discontinuities  (shocks) . 


The  same  techniques  used  to  solve  the  equations  (18)  for  a 
Cartesian  grid  can  now  be  applied  to  the  body  conforming  grid.  It  is 

important  to  remember  that  F  and  G  contain  derivatives  in  the  viscous 
and  heat  transfer  terms  that  must  be  properly  evaluated  in  the  trans¬ 
formed  system. 

Considerable  effort  went  into  defining  the  governing  equations  in 

terms  of  x,,  x  .  ye  and  y_  rather  that  £  ,  n  ,  £  and  n  .  This  was  done 
i  v  i  V  X  X*  *y  'y 

because  typically  the  numerical  approximations, 

£  =A£| 

*x  Ax  I Ay=0 
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lor  example,  are  not  easily  determined  from  a  body  conforming  computa¬ 
tional  grid.  Whereas,  the  approximation, 

Ax  I 

X(  =  A{l  A»7=0 

is  easily  determined  from  the  known  computational  grid. 

Finite  Difference  Approximations 


Finite  difference  solution  procedures  require  the  substitution  of  a 
finite  difference  approximation  in  place  of  the  differential  expression. 
The  solution  proceeds  across  the  region  of  interest  by  taking  small 
steps  in  the  spatial  direction  at  small  steps  in  time.  The  sise  of  the 
spatial  and  temporal  steps  are  critical  both  for  accuracy  of  the  solu¬ 
tion  and  stability  of  the  finite  difference  approximation.  More 
information  on  both  the  grid  (spatial  steps)  and  stability  will  be 
provided  later. 

To  see  how  the  finite  difference  approximations  are  determined, 
consider  the  definition  of  the  partial  derivative  for  f(x,y)  at  x=xo, 

y=y0» 

8f(*.y)  _  lia  f (xo+8x>y0)  -  f(xo'yQ) 

8x  Ax-*0  - - , - 

Ax 

Therefore,  in  the  limit  as  Ax  approaches  zero  the  approximation  ap¬ 
proaches  the  exact  value.  A  more  enlightening  method  of  developing 
finite  difference  approximations  is  derived  from  using  the  Taylor-series 
expansion  for  f(x+Ax,y): 

f  (x+Ax,y)  =  f  (x.y)  +  Ax  +  -|f 

9x 


t  a3f  Ax3 
8x3  31 

Solving  for  8f/8x  provides: 


(25) 


di  _  f (x+Ax ,y)  -  f(x,y)  8  f  Ax  + 

8x  Ax  n  2  2 !  +  ' 

Ox 


(26) 


The  truncation  error  associated  with  approximating  the  partial  deriva¬ 
tive  by  the  first  term  on  the  right  of  the  equal  sign  is  the  sum  of  the 
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remaining  terms  on  the  right-hand  side.  This  is  described  as  a  finite 
difference  approximation  to  the  partial  derivative  that  is  of  the  order 
of  (0)  Ax,  or  'first-order  accurate.*  Equation  (26)  is  also  referred  to 
as  a  'forward*  difference.  A  'backward*  difference  is  obtained  by  ex¬ 
panding  f(x-Ax,y)  in  a  Taylor-series  such  that; 

f (x-Ax,y)  =  f  (x,y)  -  Ax  +  ~  - 

ox  ** 


then 


03  f  Ax3 


0x 


3  3! 


(27) 


j»  .  (28) 

0x 

which  is  again  0(Ax)  or  first-order  accurate.  A  higher  order  accurate 
approximation  may  be  obtained  by  subtracting  equation  (27)  from  equation 
(25).  This  provides; 


3  3 

f  (x+Ax,y)-f (x-Ax,y)  =  2^  Ax  +  2®  f  Ax 


0x 


3  3l  +  *  *  * 


Solving  equation  (29)  for  the  partial  derivative  provides; 

8f  =  ifr-tx.,)  -  4  0(ix)2 

0x  2  Ax  v  ' 


(29) 


(30) 


Equation  (30)  is  a  'central*  difference  and  is  'second-order  accurate." 
Note  that  the  central  difference  approximation  does  not  use  the  value  of 
the  variable  at  the  point  (x,y)  of  interest. 

Finite  differences  for  second  order  (and  higher)  derivatives  can 
also  be  determined.  If  equations  (25)  and  (27)  are  added  together  the 
result  is; 

f (x+Ax,y)  +  f (x-Ax,y)  =  2f(x,y)  + 


<4*)2  + 

ax’  2! 


0(Ax) 


Solving  for  the  second  derivative; 

82*  _  f(xt-Ax.y)  -  2f(x,y)  ♦  f(x-Ax,y)  + 


0x 


(Ax)* 


0(Ax)* 


(31) 
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Equation  (31)  is  a  second-order  accurate  central  difference  for  the 
second  derivative.  These  are  but  a  few  of  the  many  possible  finite  dif¬ 
ference  approximations.  A  good  summary  of  forward,  central  and  backward 
differences  for  first-through  fourth-order  differences  is  provided  in 
Kef.  6.  Similar  finite  difference  approximations  exist  for  derivatives 
with  respect  to  y  and/or  time.  Mixed  partial  derivatives  are  also  pos¬ 
sible  and  a  summary  table  is  provided  in  Ref.  5.  The  spatial  step  sixes 
are  determined  by  the  computational  grid,  and  the  temporal  step  sise  is 
most  often  determined  by  stability  considerations. 

The  archetypal  partial  differential  equation  for  demonstrating 
finite  difference  methods  is  the  one-dimensional  linear  hyperbolic  con¬ 
vection  equation; 

+  c  =  0  c  >  0  is  constant 

Substituting  a  forward  difference  approximation  for  time  and  a  central 
difference  for  space  provides; 

u(t+At,x)  -  u(t,x)  u(t,x+Ax)  -  u(t,x-Ax) 

At  +  C  2  Ax  =  U* 

Solving  for  u(t+At,x); 

u(t+At,x)  =  (c/2) (j^) [u(t,x-Ax)  -  u(t,x+Ax)]  +  u(t,x)  (32) 

Equation  (32)  is  called  an  explicit  method  (Euler  explicit  to  be  exact) 
because  the  value  at  the  next  time  step  is  explicitly  known  in  terms  of 
values  at  the  current  time  step.  The  solution  would  be  obtained  by 
dividing  the  domain  into  small  spatial  increments  (grid) ,  applying  the 
appropriate  boundary  conditions  and  iterating  in  time  with  small  tem¬ 
poral  increments.  Explicit  schemes  typically  have  very  restrictive  time 
step  requirements  in  order  to  avoid  small  errors  from  amplifying  and 
destroying  the  solution.  In  fact,  a  von  Neumann  stability 

analysis^ ’7 indicates  that  the  explicit  Euler  method  is  unstable  for  any 
time  step. 

If  equation  (32)  were  rewritten  with  the  spatial  difference 
evaluated  at  the  next  time  step,  rather  than  the  current  time  step,  we 
obtain; 

u(t+At,x)  =  (c/2) (j^) [u(t+At,x-Ax)  - 
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u(t+At,x+Ax)]  +  u(t,x)  (33) 

Equation  (33)  is  called  an  Euler  implicit  method.  A  stability  analysis 
indicates  the  method  would  be  stable  for  any  time  step;  however,  a  set 
of  simultaneous  equations  must  be  solved  at  each  iteration.  Implicit 
schemes  typically  are  very  stable  and  allow  relatively  large  time  steps 
that  are  limited  by  truncation  errors  rather  than  stability. 
Unfortunately,  very  large  computational  domains  require  the  solution  of 
a  large  set  of  simultaneous  equations  which  may  take  considerable  com¬ 
puter  time. 

It  is  obvious  that  the  large  number  of  possible  finite  difference 
approximations  available  for  a  given  partial  differential  equation 
provides  for  an  almost  endless  number  of  possible  solution  schemes. 

Many  such  methods  are  described  in  Refs.  5  and  7.  New  methods  and 
variations  of  the  old  methods  will  undoubtedly  continue  to  evolve  for 
some  time  to  come.  Any  new  or  modified  solution  method  must  be  accom¬ 
panied  by  a  stability  analysis  to  determine  the  allowable  step  size. 
Unfortunately  the  mathematics  associated  with  such  an  analysis  is 
prohibitively  difficult  for  complex  partial  differential  equations. 

Often  an  approximate  stability  criterion  is  determined  from  analysis  of 
a  simplified  version  of  the  governing  equation,  and  then  the  final 
stability  criterion  is  determined  by  trial-and-error. 

The  solution  procedure  chosen  for  this  study  is  the  MacCormack  ex¬ 
plicit  method,  described  in  Ref.  8.  The  advantages  associated  with  this 
method  are  that  it  is  relatively  simple,  is  easily  vectorized  on  super¬ 
computers  and  provides  a  good  solution  at  flow  discontinuities  (shocks) . 
The  primary  disadvantage  is  that  the  restriction  on  the  allowable  time 
step  requires  many  iterations  for  a  converged  solution. 

MacCormack’s  method  is  called  a  two-step  or  a  predictor-corrector 
procedure.  The  first  step,  or  predictor,  is  a  first-order  forward  dif¬ 
ference  in  both  time  and  space.  The  second  step  (or  corrector)  is  a 

8 

first  order  backward  difference  in  space.  The  net  result  is  said  to  be 
second-order  accurate  in  both  time  and  space.  The  exact  derivation  of 

MacCormack’s  method  is  obscure.  Hankey  provides  the  following  finite 
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difference  "MacCormack*  representation  for  the  linear  convection  equa¬ 
tion; 

Predictor 

uCfc;x)  +  cf«0bx+Ax)  r  u(.t>x))=0  f34> 

At  v  Ax  '  v  J 


Corrector 


u(t+at,x)  -  u( 
At 


+  c(j(ti«l  -..“C ^ix.-&x.}) 


(c/2)  C“ 


[)  =  o. 


Where  u  represents  an  intermediate  (predicted)  solution.  An  inter¬ 
mediate  value  of  the  flow  variable  is  determined  at  each  point  in  the 
grid  by  solving  the  predictor  equation; 

u(t+At,x)  =  u(t,x)  -c(At/Ax) [u(t,x+Ax)-u(t,x)] .  (36) 

Solving  the  corrector  equation  (35)  for  u(t+At)  by  making  use  of  the 

definition  of  u  from  the  predictor  equation  (34)  provides; 

u(t+At,x)=(l/2){u(t,x)  +  u(t,x)  - 

(At/Ax) [u(t,x)  -  u(t,x-Ax)]>.  (37) 

Equations  (36)  and  (37)  are  the  accepted  predictor  and  corrector  expres¬ 
sion  for  the  MacCormack  method. 

The  corrector  equation  represents  a  forward-time,  backward -space 
difference.  However,  it  is  not  clear  where  the  last  term  in  equation 
(35)  comes  from.  Ref.  5  indicates  that  MacCormack ’s  scheme  is  a  varia¬ 
tion  of  the  two-step  Lax-Wendroff  procedure,  which  uses  averaging  and 
half-step  techniques.  Therefore,  a  possible  finite  difference  inter¬ 
pretation  for  the  corrector  step  that  would  also  reduce  to  equation  (37) 


(1/2)  At 


.x)  -  ul 


L]  =  0 


k 
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This  provides  a  true  backward  difference  for  the  spatial  derivative  (of 
the  predicted  variable)  but  assumes  that  the  average  of  the  original 
variable  and  the  predicted  variable  constitutes  half  of  the  time  step 
between  the  old  and  new  iteration. 

Recall  that  the  vector  form  of  the  Navier-Stokes  equations  (18) 


was; 

dU  dP  8G  _  n 
8t  +  0x  +  8y  "  U- 

In  order  to  simplify  the  notation  the  following  convention  will  be  fol¬ 
lowed  in  the  preceding  equations : 

n  represents  the  current  time  step, 
i  represents  the  current  x  location, 
j  represents  the  current  y  location, 

At  is  the  time  increment  between  n  and  n+1, 

Ax  is  the  spatial  increment  between  i  and  i+1  or 
i  and  i— 1 , 

Ay  is  the  spatial  increment  between  j  and  j+1  or 
j  and  j-1. 


Therefore,  the  forward  time,  forward  space  representation  for  the 
governing  equation  in  MacCormack’s  predictor-corrector  form  is: 


=  Un  . 
i.J  1.3 


(At/Ax) (F 


i+1,  j 


F?  .)  - 

i.J 
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(At/Ay) (Cjfj+i 


-  Gn 


i.J 


)• 


(38) 


Equation  (38)  is  used  to  determine  an  intermediate  (predictor)  value  for 
the  flow  variables  at  each  point  in  the  computational  grid.  The  correc¬ 
tor  step  is  written: 

u "+*  =  (1/2)  [U?  .  +  (At/Ax) (f°+{-  PD+1  ) 

x  I  J  A  I  J  A  I  J  A  >  J  A  A  >  J 

-  (Jt/4x)(8n41-  6”*1  ,)] 

where  P  and  ft  are  computed  using  the  intermediate,  or  predictor,  values 
of  the  vector  ft. 

MacCormack  attempted  a  stability  analysis  of  a  two-dimensional 
linearized  version  of  the  governing  equations.  After  some  rationaliza¬ 
tion  and  much  mathematics  he  arrived  at  the  stability  criterion; 

At  <  [Ax/(lul  +  Ivl  +  2c) (39) 


A  one-dimensional  analysis  provided  the  Courant-Friedrichs-Lewy 
(CFL)9  condition; 

At  <  [Ax/(lul  +  c)].  (40) 

For  small  spatial  grid  increments  equation  (39)  is  much  more  restrictive 
than  equation  (40).  UacCormack  ignores  viscous  and  heat  transfer  effects 
by  stating  they  will  not  affect  stability  providing; 

(M/p) (At/Ax2)  <  1 

and 


k (At/Ax2)  <  1. 

In  practice,  these  criteria  may  become  dominant  in  the  high  gradient 
regions  of  a  boundary  layer  where  very  small  spatial  increments  are 
defined.  Ref.  5  suggests  an  empirical  formula  for  At  that  is  based  on 
the  CFL  and  grid  mesh  Reynolds  number. 

Strong  flow  discontinuities  (shocks)  can  produce  large  gradients 
across  grid  cells.  The  resulting  numerical  oscillations  may  be  large 
enough  to  cause  the  solution  procedure  to  fail.  To  help  avoid  these 
numerical  anomalies  most  solution  procedures  provide  some  method  of 
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numerical  smoothing  (also  called  damping  or  artificial  viscosity).  The 
objective  is  to  produce  the  necessary  moderation  of  the  numerical  oscil¬ 
lations  without  affecting  the  physics  of  the  flow.  Usually  smoothing 
corrections  of  the  order  of  the  truncation  error  or  smaller  are 


8 

employed.  MacCormack  has  developed  a  fourth  order  smoothing  term  that 
is  proportional  to  the  second  derivative  of  the  pressure.  The  form  of 
the  smoothing  is; 


a  At  (Ax)3 


where  a  is  a  constant  such  that  0  <  a  <  0.5. 

For  most  flow  problems  this  smoothing  term  will  be  very  small  except  in 
the  region  of  large  pressure  gradients,  such  as  near  strong  shocks. 

Some  smearing  of  the  shocks  over  two  or  more  grid  cells  is  inevitable 
with  this  approach.  Different  solution  algorithms  and/or  adaptive  grid 
routines  may  help  minimize  the  smoothing  effects. 


Turbulence  Models 


The  Navier-Stokes  equations  model  any  fluid  flow  that  meets  the 
assumptions  made  in  their  derivation.  Unfortunately,  turbulent  boundary 
layers  present  a  situation  which,  while  theoretically  predictable  via 
the  Navier-Stokes  equations,  is  not  practical  to  do  so  because  of  the 
size  and  the  random  frequencies  of  the  turbulent  structures  in  a  high 
Reynolds  number  flow.  Numerical  resolution  of  these  turbulent  struc- 

tures  would  require  such  a  fine  grid  (10  points  per  cubic  centimeter  ) 
that  it  is  said  that  the  fastest  computer  would  take  many  years  to  solve 

the  relatively  simple  problem  of  turbulent  flow  in  a  pipe1 . 

Since  direct  calculation  of  the  small  scale  turbulent  flow  struc¬ 
tures  was  not  a  viable  option,  researchers  found  ways  of  approximating 
the  mean  flow  characteristics  of  turbulent  boundary  layers.  The  result¬ 
ing  equations  are  the  time-averaged  or  Reynolds-averaged  Navier-Stokes 
equations.  The  solution  approach  is  to  solve  the  Reynolds-averaged 
equations  for  the  mean  flow  (large  scale  flow  features)  while  modeling 
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the  fluctuating  effects  (small  scale  flow  structures)  by  some  other 
method . 

In  the  classical  Reynolds  averaging  the  primary  flow  variables  are 
replaced  by  the  sum  of  a  time-averaged  term  and  a  fluctuating  term.  For 
example,  if  f  is  an  arbitrary  flow  variable,  then 

f  =  f  +  f  , 


where  I  is  the  steady  or  time-averaged  value; 


t  +At 

o 


f  dt, 


and  f  is  the  deviation  from  this  mean.  At  is  a  time  increment  of  suf¬ 
ficient  length  to  provide  a  statistically  significant  result.  By 
definition,  the  time  average  of  the  fluctuating  part  is  zero.  That  is; 

>'  ■  it  **4t «'  * .  o. 

O 

However,  the  time-average  of  the  products  of  the  fluctuating  properties 
is  not  zero  (TT  >  0) . 

To  demonstrate  the  Reynolds  averaging  process,  the  continuity  equa¬ 
tion  (10)  for  one  dimension  is  written  as; 

=  °.  (41) 


Substituting  p  =  p  +  p  and  pu  =  (p  +  p  ) (u  +  u  )  =  pu  +  p  u  +  pu  + 
p  u  into  equation  (41)  yields; 


PR  +  §£  +  $£2  +  +  &L.  +  j£JL  =  0. 

dt  9t  9x  9x  9x  9x 
Taking  the  time  average  of  equation  (42)  provides; 


(42) 


+  l£_  +  $£*  +  +  +  Pj-1 L  =  0 

9t  9t  9x  9x  8x  9x  ’ 

where  the  only  non-zero  terms  are; 


bp  bpu  bp  u 
9t  +  9x  *  9x 


=  0 


or  since  f  =  f  (see  Ref.  5  or  Ref.  10  for  a  more  complete  derivation) 
this  is  the  same  as; 
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+  «.  ®£_ 
9t  9x  dx 


=  0. 


i 


(43) 


I 


Comparing  equation  (43)  and  (41)  reveals  that  the  continuity  equation 
using  the  original  variables  differs  from  the  continuity  equation  using 


the  time-averaged  variables  only  by  the  addition  of  a  function  of  p  u 
The  Reynolds-averaged  form  of  the  momentum  equation  and  the  energy  equa-  I 

tion  follow  a  similar  procedure.  In  each  case  the  Reynolds-averaged 
equations  using  time-averaged  variables  are  equal  to  the  original  equa¬ 
tions  plus  some  terms  containing  time-averaged  products  of  fluctuating 
variables.  For  example,  the  momentum  equation  in  the  x  direction  for  I 


two-dimensional  flow  may  be  written  as 

8_ 

dx 


1 1 


a  _  '  /  &  _  _  _ _ 

^0>U  +  P  u  )  +  ^(/Juu  +  U/7  u  )  +  j^(puv  +  up  V  )  = 


9_ 
dyf 


9p  9  ,  ro9u  . _._w9u  9vN, 

"  dx  +  _  (2/3)(ftT  +  ar)]  ' 


9x^L  9x 


9x  9y' 


up  u  -  pu  u 


,9u  9v 


-  U  r  /U  U  UTv - -  -  -i 

+  371X07  +  37)  ~  u  -/»uv  -  /»uvj 


(44) 


Using  the  same  definition  for  as  for  (equation  (14))  the  right- 
hand  side  of  equation  (44)  becomes; 

0^(TXX  ~  ”  P  u  u  )  + 

0  0  *  0  0  0  0  0 
9y^xy  "  '  P  u  v  ^ 

The  time  averages  of  the  fluctuating  variables  are  said  to  look  like 


stress  terms,  -^u  v  and  -pu  u  are  specifically  identified  as 
"Reynolds  stresses."  Note  that  p  is  not  time  averaged  in  this  defini¬ 
tion  because  the  classical  definition  for  Reynolds  stresses  was  derived 


for  incompressible  flow1.0  These  apparent  stresses  provide  a  measure  of  I 

the  momentum  transport  due  to  the  turbulent  motions  in  the  boundary 
layer  Applying  the  same  averaging  techniques  to  the  energy  equation 

produces  apparent  heat-flux  terms&  (pi  u  ,  p  T  u  ,  etc.). 

I 
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The  time  averaging  process  provides  equations  for  the  large  scale 
features  of  the  flow  but  at  the  expense  of  introducing  several  addi¬ 
tional  unknowns.  Inconpressible  analysis  greatly  simplifies  the 

Reynolds  averaged  equations  because  all  terns  containing  p  would 
vanish.  Coapressible  flow  equations  are  sinplified  somewhat  by  using 

order-of -magnitude  arguments1.  ’5  In  recent  years  it  has  become  popular 
to  develop  Reynolds-averaged  equations  using  mass-weighted  averaging  for 
the  flow  and  thermal  variables.  This  approach  produces  relatively 
simple  momentum  and  energy  equations  without  making  simplifying  assump¬ 
tions.  However,  much  more  complex  expressions  are  obtained  for  the 
stresses,  and  order-of -magnitude  arguments  are  required  in  order  to 
evaluate  them.  Ref.  5  provides  a  good  comparison  between  classical 
time-averaging  and  mass-averaging. 

Determining  what  to  do  about  the  remaining  new  terms,  and  how  to  do 
it  is  referred  to  as  the  "closure  problem*  for  Reynolds-averaged  Navier- 
Stokes  solutions.  In  the  development  of  the  original  Navier-Stokes 
equations  researchers  looked  to  the  properties  of  the  fluid  for  closure 
information  (e.g.  the  perfect  gas  relationships).  The  search  for  ways 
to  evaluate  the  products  of  the  fluctuating  properties  is  referred  to  as 
"turbulence  modeling"  and  relies  heavily  on  experimental  observations. 

Hany  different  turbulence  models  with  wide  degrees  of  complexity 
have  been  proposed.  Probably  the  simplest,  and  the  only  one  to  be  dis¬ 
cussed  here,  is  based  on  the  Boussinesq1  assumption  that  the 

Reynolds  stresses  are  related  to  the  laminar  stress  terms  through  an 
"eddy"  viscosity  term  (e). 


-pn 


■4 


— — '  ro9u  /o/ow^U  ®VM 
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~Py  v  =  ~  +  0^)] 

Note  that  the  value  for  e  may  be  different  in  each  equation;  however, 
the  convention  is  to  assume  they  are  equivalent.  The  local  mean  density 
is  used  for  the  eddy  viscosity  definition.  The  other  fluctuating  terms 


22 


are  assumed  to  be  snail  and  are  neglected  in  the  definition  of  the  tur¬ 
bulent  shear  stress.  These  simplifications  appear  to  work  well  for 

attached  boundary  layers  at  Mach  numbers  up  to  5S .  The  obvious  ad¬ 
vantage  of  this  definition  is  that  the  total  viscosity  is  then  the  sum 
of  the  thermodynamic  viscosity  fl  and  the  eddy  viscosity  6. 

In  a  similar  manner  an  "eddy*  conductivity  (k.)  may  be  defined  to 

t 

approximate  the  turbulence  effects  on  the  energy  equation; 


-'V  T  =  kt  §7 

Fortunately,  experiment  has  shown  that  the  eddy  conductivity  can  be  re¬ 
lated  to  the  eddy  viscosity  through  the  turbulent  Prandtl  number 

(PrJ1;’5 


k.  =  c  e/Pr, . 
t  p  '  t 

Where  Pr^  is  assumed  to  be  a  constant  with  a  value  near  1.0.  Closure  of 

the  Reynolds  averaged  Navier-Stokes  equations  now  depends  solely  on  the 
evaluation  of  the  eddy-viscosity .  Unfortunately  the  eddy-viscosity  is 
not  a  simple  physical  property  of  the  fluid  as  is  thermodynamic  vis¬ 
cosity  fi. 


Algebraic  Eddy  Viscosity  Model 


Algebraic  or  zero-equation  (meaning  no  partial  differential 
equations)  models  for  e  invariably  incorporate  the  Prandtl  mixing  length 
assumption; 


i  UU  i 

6  =  p  *- 


du, 

wiH  . 

9y 


Where  H  is  the  transverse  distance  over  which  a  turbulent  structure 


maintains  its  original  momentum,  u  is  the  velocity  in  the  primary  flow 
direction  (x) ,  and  y  is  transverse  to  u.  The  remainder  of  this  dis¬ 
cussion  will  be  on  the  Baldwin-Lomax1  1  model  used  in  this  study.  This 

description  is  for  two-dimensional  flow  and  the  variables  u,  v,  x,  and  y 
are  as  defined  above  but  the  *  will  be  dropped. 


A 
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The  Baldwin-Lomax  model  is  a  two-layer  algebraic  model  with  an  in- 
teraittency  factor.  Two-layer  means  that  there  are  two  separate 
equations;  one  for  the  inner  layer  that  forces  the  eddy  viscosity  to 
approach  zero  at  the  wall,  and  another  equation  for  the  remainder  of  the 
boundary  layer.  The  intermittency  factor  is  an  experimentally  derived 
expression  that  accounts  for  the  decreasing  frequency  of  encountering 
turbulent  flow  structures  in  the  flow  outside  of  the  boundary  layer. 

An  eddy-viscosity  coefficient  is  defined  at  each  x  location.  The 
equation  for  the  inner  eddy-viscosity  (e^)  is; 


Cj  =  p  fc2  Iwl .  (45) 

Where 

l  =  kyD, 

k  is  the  von  Karman  constant  =0.40 
D  is  the  van  Driest  damping  function; 

D  „  1  -  .<-*>*), 


A  is  a  sublayer  constant  =  26, 

«.  fPr )‘/2 

y  =  y  (-^— '  )  I  ^  (subscript  w  refers 

to  the  solid  surface), 

.  .  . ,5v  9u. . 

*  'fe'ar5'  ■ 


The  von  Karman  constant  provides  the  proper  scaling  for  £-  near  the  solid 
surface.  The  van  Driest  function  provides  for  the  transition  from  the 
viscous  sublayer  to  the  fully  developed  turbulent  outer  region. 
Therefore,  £.  varies  from  zero  at  the  surface  to  essentially  a  linear 
function  at  large  distances  from  the  surface.  Iwl  is  the  magnitude  of 
the  vorticity  which  is  typically  very  large  near  the  surface,  then 
decreases  rapidly  through  the  boundary  layer. 

The  equation  for  the  outer  eddy  viscosity  (eo)  for  attached  flow 


over  a  solid  surface  is; 

where , 


Co-'K  Ccp  Fwake  Fkleb ‘ 


(46) 


J 


J 


« 


J 
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K  is  the  Clauser  function  =  0.0168, 

C  =  1.6, 
cp 

F  ,  =  Y  F 

wake  max  max 


(47) 


The  Clauser  constant  is  related  to  the  Reynolds  number  based  on  boundary 
layer  thickness6  and  C  is  a  constant  applied  to  this  method  to  match 


1  2 

the  Cebeci-Smith  results. 

F(y)  =  y  Iwl  D 


F  is  determined  from  the  equation; 
max 


by  locating  the  y  value  at  which  F(y)  is  a  maximum  for  a  given  x  loca¬ 
tion.  Y  is  the  value  of  y  that  produces  F  This  F(y)  function  is 

reported11  to  be  a  much  more  reliable  indicator  than  trying  to  find  the 

edge  of  the  boundary  layer  as  in  previous  (e.g.  Cebeci-Smith  )  methods. 

The  edge  of  the  boundary  layer  is  not  a  clearly  defined  location. 
The  random  encounter  with  turbulence  above  the  mean  edge  of  the  boundary 
layer  is  referred  to  as  intermittency .  Interaittency  effects  are  ac¬ 
counted  for  by  the  Klebanoff  intermittency  factor  (F^^) ; 


where  C 


kleb 


kleb 
=  0.30. 


=  [  I-  *  5-5(Cklebf- 


max 

^kleb  *s  a  ^unc^on  such  that  the  value  of  F^j^  is 


near  1  at  small  values  of  y  but  diminishes  rapidly  as  y  extends  beyond 
the  edge  of  the  boundary  layer. 

With  a  slightly  different  definition  for  F^^  the  outer  eddy  vis¬ 


cosity  formulation  is  reported1 1  to  be  valid  for  separated  flow  and 
wakes.  The  separated-f low/wake  formulation  for  Fw&j{e  is: 


where 


F _ i__  =  (C_,_  Y 


wake 


wk  max 


Udif^F 


max' 


(48) 


Cwk  is  another  empirical  constant  =  0.26, 


Udif 


Is  \  1  f  2  2  X \  1 

'max  '  'min’ 
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and  Faax  is  as  defined  above.  For  separated  flow  over  a  solid  surface 

the  second  term  for  U,..  is  sero,  and  for  wake  flow  D  =  1  in  the  F 

dif  ’  max 

equation.  To  avoid  having  to  locate  the  point  of  separation,  if  any, 

Ref.  11  suggests  calculating  both  ways  (equations  (47)  and  (48)) 

and  using  the  saallest  value. 

For  a  given  x  location,  both  an  inner  and  an  outer  eddy  viscosity 
coefficient  have  been  defined  (equations  (45)  and  (46)  respectively). 
Obviously  both  equations  will  provide  a  value  for  e  for  any  value  of  y. 
Ref.  11  suggests  calculating  both  values  and  switching  from  using  the 
inner  value  (e^)  to  the  outer  (eo)  value  at  the  smallest  value  of  y  for 

which  e.  >  e  .  Note  that  this  criterion  is  not  the  same  as  simply 
i  o 

choosing  the  smaller  of  the  two  values. 

It  is  important  to  remember  that  the  eddy  viscosity  model  coeffi¬ 
cients  are  based  on  experimental  observations,  usually  on  flat  plates 
and  at  subsonic  Mach  numbers.  The  values  for  the  constants  provided  are 

all  subject  to  debate1, 3  ’ 1 4  ’ 1 5  particularly  in  regions  of  high  pressure 
gradients,  wakes,  and  at  high  transonic  or  supersonic  Mach  numbers. 

Also,  many  researchers  are  very  critical  of  using  the  algebraic  equation 

16  17  18 

model  for  flow  conditions  that  are  not  fully  attached.  *  ’ 

The  reason  for  emphasis  on  the  Baldwin-Lomax  procedure  and  its  in¬ 
herent  deficiencies  will  become  apparent  in  the  next  section  where  a 
modification  of  the  Baldwin-Lomax  approach  that  allows  solutions  for 
flow  with  massive  shock-induced  separation  is  described.  First, 
however,  a  discussion  on  boundary  conditions  and  computational  grids  is 
necessary. 

Boundary  Conditions 

A  discussion  on  boundary  conditions  should  probably  begin  with  a 
discussion  on  existence,  uniqueness  and  *well-posedness . *  Ref.  15 
points  out  that  there  is  no  rigorous  mathematical  theorem  to  define  the 
criteria  to  ensure  existence  and  uniqueness  of  solutions  for  the  Navier- 
Stokes  equations.  Ref.  7  points  out  that  uniqueness  is  particularly 
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uncertain,  and  refers  to  stall  hysteresis  and  other  examples  where  dif¬ 
ferent  flow  conditions  may  exist  for  identical  boundary  conditions. 

7  19  2  0 

Existence  of  solutions  appears  to  be  somewhat  less  of  a  problem. ’  ’ 

A  problem  is  said  to  be  well-posed  if  any  small  change  in  the 
boundary/initial  conditions  produces  a  corresponding  small  change  in  the 
1  9 

solution.  Kef.  20  suggests  that,  for  the  one-dimensional,  incompres¬ 
sible  form  of  the  Navier-Stokes  equations,  well-posed  conditions  are  (1) 
known  fixed  initial  conditions  for  each  variable,  and  (2)  Dirichlet 
boundary  conditions.  Dirichlet  boundary  conditions  are  one  of  the  three 
classes  of  boundary  conditions  defined  in  Ref.  2  as: 

Dirichlet  -  Values  for  the  variables  are 
specified  on  the  boundary. 

Neumann  -  The  normal  derivative  of  the  variable 
is  specified  on  the  boundary. 

Robins  (or  mixed)  -  A  linear  combination  of 

Dirichlet  and  Neumann  boundary  conditions. 

Solid  surface  boundary  conditions  for  the  momentum  equation  vari¬ 
ables  (u,  v)  are  exceedingly  simple,*  velocity  relative  to  the  surface 
equals  zero.  The  energy  equation  variable  (e)  is  usually  determined  by 
computing  the  temperature  at  the  surface,  then 
e  =  C  T 

W  V  w 

where  the  subscript  w  refers  to  the  wall  or  surface.  Remember  that  the 
velocity  at  the  wall  is  zero.  The  two  most  common  surface  conditions 
for  temperature  are: 

T  =  constant 

w 


<8>.  -  °- 

Classical  heat  transfer  equations  could  be  employed  for  more  difficult 
cases.  Equilibrium  temperature  is  reached  very  quickly  for  metallic 
models  in  continuous  flow  wind  tunnels.  Therefore,  the  solution  proce¬ 
dure  discussed  in  this  report  uses  a  turbulent  flow  recovery  factor  to 

estimate  wall  temperature1.0 
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The  continuity  equation  variable  (/>)  remains  to  be  defined  at  the 
solid  surface.  This  is  usually  accomplished  by  applying  the  experimen¬ 
tally  observed  characteristic  of  the  pressure  gradient  through  the 
boundary  layer  perpendicular  to  the  surface; 

<§?>.  -  <«> 
(Ref.  2  obtains  this  result  by  examining  the  momentum  equation  normal  to 
the  surface.)  Using  the  equation  of  state  in  equation  (40)  provides  a 
finite  difference  approximation  for  density  at  the  wall; 

T2 

P*  ~  pl  ^T  ^ ' 

w 

The  subscript  2  refers  to  the  first  grid  level  above  the  surface.  Note 
that  p is  a  compatibility  condition,  not  a  true  boundary  condition  and 

the  smaller  the  first  grid  step,  the  better  this  approximation. 

A  *C*  type  computational  grid  (Figure  1)  was  used  for  the  solutions 
in  this  study  (grid  generation  is  discussed  in  the  next  part) .  As  a 
result,  there  is  a  branch  cut  extending  from  the  trailing  edge  of  the 
airfoil  to  the  downstream  boundary.  The  upper  and  lower  parts  of  the 
grid  have  common  points  along  the  cut  line.  Values  of  the  flow  vari¬ 
ables  along  the  cut  line  are  defined  as  the  average  of  the  value  at  the 
first  level  above  and  the  first  level  below  the  cut  line. 

The  most  difficult  and  controversial  boundary  conditions  arise  at 
the  "far-f ield"  boundaries.  Problems  occur  because  the  computational 
domain  for  the  numerical  solution  is  necessarily  limited  in  both  space 
and  time.  For  example,  an  F-15  flying  at  30,000  feet  has  a  far-field 
boundary  (the  earth)  about  500  fuselage  lengths  away.  The  most  generous 
numerical  boundary  would  have  to  be  an  order  of  magnitude  or  more 
smaller.  The  task  then  is  to  devise  physically  accurate,  stable  and, 
hopefully,  uncomplicated  boundary  conditions  for  the  flow  variables  at 
these  artificial  boundaries. 

In  supersonic  flow  disturbances  cannot  propagate  upstream  from  the 
body,  therefore,  boundaries  with  flow  into  the  computational  domain  must 


28 


be  at  free  stream  conditions  (e.g.,  p  =  pm,etc.) .  Likewise,  distur¬ 
bances  downstream  of  the  body  cannot  propagate  upstream  to  influence  the 
flow  over  the  body.  As  a  result,  the  simple  procedure  of  extrapolating 
the  values  of  the  variables  from  the  interior  of  the  domain  to  the 
downstream  boundary  provides  satisfactory  results. 

Subsonic  flow  conditions  present  the  most  difficulties  for  estimat¬ 
ing  the  proper  values  for  the  flow  variables  at  the  far  field 

computational  boundaries.  Recent  literature  focuses  on  three  primary 
approaches.  The  most  complex  of  these  involves  defining  the  primary 
flow  variables  in  terms  of  characteristic  variables  derived  from  two- 

.21 

dimensional,  inviscid  equations.  A  form  of  this  approach,  referred  to 

as  Riemann  invariant,  is  used  to  define  the  upstream  boundary  condi¬ 
tions  in  the  code  used  for  this  investigation.  Table  1  provides  a 
summary  of  the  Riemann  invariant  procedure.  Next  in  level  of  complexity 

.23 

is  the  non-reflecting  boundary  conditions.  As  the  name  implies,  the 
objective  is  to  allow  disturbances  produced  in  the  solution  domain  to 
"pass  through"  the  boundary  and  not  be  reflected  back  to  contaminate  the 
solution.  This  procedure  was  not  used  in  the  current  study,  and  the 
interested  reader  is  referred  to  Ref.  23.  The  simplest  boundary  condi¬ 
tions  use  free  stream  values  for  the  primary  variables  at  the  upstream 
boundary,  and  "no-change"  conditions  at  the  downstream  boundary.  No¬ 
change  simply  means  that  the  value  of  the  variable  at  the  boundary 
equals  the  value  at  the  first  computed  grid  point  inside  the  boundary. 
This  approach  over -specif ies  the  boundary  conditions,  but  produces 
2  4 

reasonable  results  when  the  boundary  is  located  sufficiently  far  from 
the  body.  Ref.  24  provides  a  good  summary  of  the  three  approaches  dis¬ 
cussed  above  and  compares  the  results  of  applying  each  for  the  solution 
of  flow  around  a  cylinder. 

With  the  numerical  boundaries  at  sufficiently  large  distances  from 
the  body  (say  50  body  lengths),  the  effects  of  the  type  of  boundary  con¬ 
dition  applied  would  be  expected  to  vanish.  However,  computer  memory, 
run  times  and  the  error  introduced  by  highly  stretched  grids  put  severe 
constraints  on  practical  locations  for  the  far  field  boundaries.  The 
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Riemann  invariant  approach  would  seem  to  provide  a  more  physically  ac¬ 
curate  representation.  Ref.  25  claims  very  good  results  with  such  an 

approach  for  a  finite  volume  solver.  However  Shang  observed  little  or 
no  difference  between  the  characteristic  method,  the  non-reflecting 
method  or  the  simple  free  stream  and  no-change  approach  when  used  with 
an  explicit  finite  difference  solver  on  flow  over  a  cylinder.  This 
author  has  used  both  the  Riemann  invariant  and  the  free  stream/no-change 
boundary  conditions  for  flow  over  an  airfoil  with  no  observable  dif¬ 
ferences.  Both  this  study  and  Ref.  24  obtained  viscous  flow  solutions 
with  far-field  boundaries  at  about  10  reference  body  lengths,  and  both 
used  a  similar  explicit  finite  difference  code. 

The  selection  of  far  field  locations  and  the  boundary  conditions  to 
be  applied  are  still  an  inexact,  controversial,  and  continuously  evolv¬ 
ing  art.  The  results  provided  in  this  report  used  the  Reimann  invariant 
boundary  conditions  for  inflow  (upstream)  boundaries,  and  the  no-change 
conditions  for  outflow  (downstream)  boundaries.  These  boundary  condi¬ 
tions  were  chosen  because  they  seem  more  physically  proper  and  because 
other  government  and  industry  personnel  conducting  similar  studies  have 

2  6 

had  good  success  (private  conversations) .  Visbal  obtained  good 
results  for  transonic  flow  over  airfoils  with  free  stream  far-field 
boundary  conditions  when  the  boundary  was  sufficiently  far  away.  The 
far-field  boundaries  for  results  presented  in  this  study  were  located  at 
approximately  15  chord  lengths  away  from  the  airfoil. 

Computational  Grid 

Finite  difference  solutions  of  the  partial  differential  equations 
governing  fluid  flow  about  arbitrary  bodies  require  the  development  of  a 
computational  grid  around  the  body.  A  grid  system  that  conforms  to  the 
shape  of  the  body  at  the  surface  has  the  very  important  advantages  of 
simplifying  the  partial  differential  equation  solution  technique  and 
simplifying  the  satisfaction  of  the  surface  boundary  conditions.  The 
character  of  the  entire  flow  field  is  determined  in  the  high  gradient 
regions  near  the  body.  Therefore,  the  development  of  a  body-conforming 
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coordinate  system  and  computational  grid  is  a  very  important  first  step 
in  the  numerical  solution  of  fluid  flow  problems. 

Algebraic  grid  generation  techniques  have  been  used  for  some  ap- 

2  7 

plications  and  have  the  advantage  of  being  relatively  simple  and  fast. 
However,  the  most  popular  techniques  use  the  solution  of  partial  dif¬ 
ferential  equations  (PDE) .  Procedures  exist  that  solve  elliptic, 
hyperbolic  and  parabolic  PDEs.  Elliptic  routines  are  most  often  used 
and  probably  are  more  generally  applicable  to  a  wider  range  of  problems. 
Hyperbolic  routines  are  relatively  fast  and  simple,  but  are  not  easily 
adapted  to  physical  constraints  in  the  computational  domain  (such  as 

internal  flows  or  multiple  bodies) .  Parabolic  routines  have  been 

recently  developed  and  show  some  promise.  The  PDE  methods  attempt  to 
produce  nearly  orthogonal  grids  with  assurance  of  no  crossover  of  ad¬ 
jacent  grid  lines. 

A  hyperbolic  grid  generation  routine  was  chosen  for  this  study  be¬ 
cause  it  is  quite  fast,  provides  nearly  orthogonal  grids  and  has  good 
user  control  of  the  grid  spacings.  The  routine  appears  to  work  well  for 
generating  a  grid  around  two-dimensional  bodies  when  there  are  no  other 
physical  constraints  within  the  domain  of  the  grid.  The  actual  grid 
generation  starts  at  the  body  at  user  defined  locations  and  marches  out 
to  an  outer  boundary. 

The  core  of  the  program  was  adapted  from  the  0-grid  procedure 
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described  by  Steger  and  Chaussee.  That  procedure  was  modified  for  this 
study  by  including  additional  dissipation  terms  and  by  changing  the  form 
of  the  dissipation  described.  The  modifications  necessary  to  produce  a 
C-type  grid  are  also  described  and  an  example  is  provided. 

It  is  desirable  that  the  distribution  of  points  around  the  body  be 
2  7. 

second-order  smooth.  Since  most  bodies  are  defined  with  a  somewhat 
random  distribution  of  points,  a  routine  was  added  to  redefine  the  coor- 

3  0 

dinates  in  a  smooth  distribution.  This  routine  also  allows  selective 
clustering  of  grid  points  in  regions  of  high  gradients. 
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In  the  grid  generation  algorithm  developed  by  Steger  and 
2  8 

Chaussee,  definitions  of  grid  orthogonality  and  the  transformation 
Jacobian  were  chosen  to  devise  a  scheme  mapping  x,y  into  £,17: 

xiKv  +  ytyv  =  0  (5°) 

Y,  -  y{  = 1'3  <51> 

Note  that  by  the  definition  of  the  Jacobian,  cell  areas  in  physical  and 
computational  domains  can  be  written; 
dA  =  dxdy  =  (1/J)  d£dtj 

Numerically  the  grid  spacings  in  the  transformed  plane,  A£  =  Aty  =  1,  so 
the  inverse  Jacobian  approximates  the  physical  cell  area.  Equations 

(50)  and  (51)  are  locally  linearized  using; 

♦ 

x  =  x  +  Ax 
♦ 

y  =  y  +  Ay 

•  ♦ 

where  x  ,y  is  a  nearby  location.  The  linearized  set  of  equations  be¬ 
come  (after  some  algebraic  manipulation) 


•  ♦ 

* 

X 

♦  ♦ 

X 

0 

V  V 

t  t 

*+* 

= 

♦  • 

y  -x 

V  V 

y 

£ 

-ji  \ 

y 

V 

A+A* 

or,  substituting  symbols  for  the  matrices  and  vectors; 

C  ^  +  B  ft  =  1  (52) 

The  specific  numerical  procedure  for  solving  equation  (52)  depends 
on  the  class  of  equation  (i.e.,  elliptic,  parabolic,  or  hyperbolic). 

The  class  can  be  determined  by  inspecting  B  'c. 


B  1 C  =  1/J 


•  •  ♦  •  ♦  *  •  ♦ 

X,  X  -  Y,  Y  X.  Y  +  Y*  X 

£  v  £  v  tv  tv 


X  Y.  +  Y  Xf 
V  t  V  t 


Ye  Y  -  X  X- 
tV  V  t 


where , 
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2 


.  .2 

J  =  +  Yc 


,-i 


Since  B  C  is  symmetric,  it  has  real  eigenvalues,  specifically; 


•  j  ♦  *  •  •  0  J  ♦  »  •  , 

X  =  *  1/J  \/(X.  X  -  Y.  Y  )  +  (X,  Y  +  Yf  X  ) 

l  ,2  '  i  V  (  7  v  V  (  r/J 

This  indicates  that  the  system  is  hyperbolic  and  can  be  marched  in  the  17 

direction  in  a  time-like  manner.  The  finite  difference  solution  scheme 

chosen  to  approximate  equation  (52)  is  centrally  differenced  in  £  and 

backward  differenced  in  17.  The  scheme  can  be  written  as; 

[I  +  B-1C  «.]  L  ,  =  L  +  B-1  I.  .  , 

(  i,j+l  i,j  i,j+l 

+  a  (V.  A.)2  R.  . 
v  1  V  i,j 

where  6 ^  represents  a  central  difference  in  the  £  direction,  i  repre¬ 
sents  the  index  in  the  £  direction,  and  j  represents  the  index  in  the  17 
direction.  assumed  to  be  known  at  the  j  level  by  a  process 

shown  below,  and  o  (V.  A.)  R.  .  is  an  added  fourth  order  dissipation 

1  1  1 ,  j 
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term  in  £  to  aid  in  stability. 

Typically  and  y’  are  approxiaatad  as; 

-  <Vi,j  - 

y£  =  (yi*l,j  -  yi-l,i)/2 

•  ♦ 

While  x  and  y  are  extracted  from  the  nonlinear  differential  equations 
7  7 

(50)  and  (51) . 


xfj  ~  (-y^  A  )/J  , 


y7  =  A  ’ 


This  allows  the  nonlinearity  to  be  maintained. 

As  mentioned  before,  the  cell  areas  at  the  j  and  j+1  levels  are 
assumed  to  be  known.  This  can  be  accomplished  in  any  number  of  ways. 


2  8 

Steger  and  Chaussee  construct  polar  grids  about  two  individual  circles 
whose  circumference  is  the  arc  length  of  the  body  to  be  grided.  The 
grids  on  these  two  circles  have  the  same  spacing  in  the  radial  direction 
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but  differ  in  grid  spacing  in  the  circumferential  direction.  One  circle 
has  points  placed  in  equally  spaced  increments  around  the  circumference 
and  the  other  circle  has  the  same  spacing  as  the  point  distribution  over 
the  body.  Cell  areas  are  extracted  from  these  polar  grids  such  that 
when  near  the  body,  cell  areas  are  calculated  using  the  non-uniform  dis¬ 
tribution  of  points,  and  when  far  from  the  body,  areas  from  the  uniform 
distributed  circle  are  used.  A  smooth  function  transitions  from  the 
inner  areas  to  the  far-field  areas  at  a  user  defined  rate. 

The  solution  algorithm  proposed  for  these  equations  exhibited  dif¬ 
ficulties  for  geometries  with  slope  discontinuities  and  regions  of 
reverse  (concave)  curvature .  Discontinuities  propagated  into  the  grid 
interior  with  undesirable  results.  Drawing  from  experience  in  other 
hyperbolic  systems  it  was  postulated  that  these  problems  could  be  cir¬ 
cumvented  without  drastically  compromising  orthogonality  by  carefully 
including  other  forms  of  dissipation.  Using  the  analogy  between  march¬ 
ing  in  time  and  marching  in  r] ,  a  time-like  dissipation  was  added  to 
insure  smoothness.  To  accomplish  this  task,  a  more  general  class  of 
integration  in  rj  was  chosen: 


A>7 


— —  ^  ^  9*7  j  +  °  fty  j+1 


When  a  =  1,  the  original  backward  differenced  scheme  is  obtained.  The 
numerical  error  term  for  this  scheme  is 


0  R 

(1/2)  (l-2a)  hrj  — —  +  higher  order  terms. 

W 

Note  that  for  a  =  1 /2  the  integration  is  a  trapezoidal  type  and  is  for¬ 
mally  second-order  in  ij.  For  a  >  l  /2  the  error  term  has  a  dissipative 
effect  in  the  r\  direction.  Rewriting  the  algorithm  in  delta  form,  such 


that  increments  in  R  (A&  =  iL  +  ^  ~  ^j)  nay  be  determined,  the  algorithm 
becomes : 


[I  +  a  B"!C  $,]($.  .  -  5.)  =  BT1  (aA°.  1+  (l-a)A.) 

s j+i  j  j  j+l  j 

a(7.A.)2l. 

v  i  i'  j 


+ 
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Adding  second-order  implicit  dissipation  to  this  algorithm  serves 
to  augment  the  amount  of  explicit  dissipation  that  may  be  added  before 
degrading  the  accuracy  of  the  method.  Therefore,  the  final  form  of  the 
algorithm  with  both  implicit  and  explicit  dissipation,  and  the  "alpha* 
parameter,  a  is: 

[I  +  a  (V.A.)  +  a  B1  C  fij (S.  ,  -  J.)  = 

L  1  v  1  1'  v  j+1  y 

B-1  (a  A0.  ,  +  (l-a)A.)  *  an  (V.A.)2  ft. . 

v  '  j'  2vll/  j 

a j  is  the  implicit  dissipation  coefficient  and  a is  the  explicit  dis¬ 
sipation  coefficient. 

The  marching  algorithm  allowed  by  the  hyperbolic  form  of  the  equa¬ 
tions  provides  an  efficient  means  for  constructing  a  computational  grid. 
A  grid  solution  begins  by  initializing  the  variables  for  grid  stretch¬ 
ing,  dissipation  coefficients,  a,  and  reading  in  the  user  defined  body 
coordinates.  Table  2  provides  an  example  cf  the  values  of  the 
parameters  used  to  produce  the  computational  grids  in  this  study.  The 
marching  procedure  in  the  r\  direction  (increasing  j)  consisted  of  four 
steps.  First,  the  areas  at  the  j+1  level  were  approximated  by  the 
method  described  above.  Then  a  system  of  equations  was  set  up  to  deter¬ 
mine  the  increments  in  x  and  y.  This  resulted  in  a  2x2  block 
tridiagonal  matrix  to  be  solved.  The  increments  in  x  and  y  were  then 
known  at  the  j  level  for  each  £  (i)  station,  and  grid  points  at  the  j+1 
level  were  computed  using  these  increments.  The  process  was  repeated 
until  the  maximum  grid  level  was  reached  at  the  outer  boundary. 

0-type  computational  grids  are  widely  used  for  a  variety  of  solu¬ 
tion  procedures  (Euler  and  full  potential  for  example).  For  viscous 
flow  applications,  however,  a  C-type  grid  is  preferred.  The  viscous 
flow  characteristics  around  an  airfoil  and  the  resulting  wake  area  be¬ 
hind  the  airfoil  require  fine  grid  spacings  close  to  the  body  and  in  the 
wake. 

The  most  obvious  change  required  in  converting  an  0-grid  routine  to 
a  C-grid  routine  is  the  establishment  of  a  downstream  boundary  spacing 
procedure  (boundary  condition) .  A  separate  subroutine  was  created  to 
compute  the  downstream  boundary  spacings.  x  is  fixed  at  the  maximum 


(downstreaa)  value  and  y  is  incremented  as  *  &i,  the  step  size  in  the 
normal  (i j)  direction.  The  terms  that  defined  the  periodic  nature  of  the 
0-type  grid  were  eliminated  from  the  tridiagonal  matrix.  Figure  1  shows 
a  comparison  of  an  0-type  and  a  C-type  grid  around  the  NLR  7301  airfoil 
analyzed  in  this  study . 

The  hyperbolic  grid  generating  procedure  is  very  fast  (a  few 
seconds  on  a  PRIME  750)  and  provides  a  grid  pattern  with  a  lot  of  user 
control.  The  routine  is  tolerant  to  variations  of  all  input  parameters 
(robust)  except  for  the  parameter  that  controls  the  transition  to  equal 
cell  areas.  Rapid  changes  in  grid  spacing  around  the  body  or  concave 
regions  may  cause  a  shock-like  discontinuity  to  develop  and  propagate 
through  the  grid.  Tolerance  to  these  conditions  may  be  improved  (at 
some  sacrifice  in  orthogonality)  by  changing  0^  ,  o^,  or  a  and/or  by  ad¬ 
ding  or  deleting  grid  points  along  the  body. 
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II.  APPLICATION  TO  A  SUPERCRITICAL  AIRFOIL 


Current  State  of  the  Art 

The  objective  of  this  Navier-Stokes  work  was  to  compute  the  two- 
dimensional,  steady-state  flow  over  a  16.5%  thick  supercritical  airfoil 
(NLR  7301)  at  a  Mach  number  of  0.807  and  an  angle  of  attack  of  0.36  de¬ 
grees.  Reynolds  number  based  on  the  chord  length  ranged  from  3  to  12  x 

10  and  boundary  layer  transition  from  laminar  to  turbulent  occurred 
naturally  at  approximately  10  percent  chord.  At  these  conditions  there 
is  a  strong  shock  on  both  the  upper  and  lower  surfaces  with  separated 
flow  extending  from  the  shocks  to  the  trailing  edge  (Figure  2) .  These 
conditions  present  a  severe  test  case  for  computational  methods. 
Successful  solutions  (solutions  that  accurately  match  experimental  data) 
for  flow  conditions  with  strong  shocks  and  extensive  separation  are 
rare.  Poor  agreement  between  experimental  and  computational  results  may 
be  caused  by  such  things  as  inadequate  computational  grids,  improper 
boundary  conditions  (which  includes  wind  tunnel  wall  effects)  and  in¬ 
valid  turbulence  models. 

Computational  grids  with  as  few  as  5,400  points  (135  x  40)  to  as 
many  as  50,000  points  (500  x  100)  were  used  in  this  study.  Even  a  two- 
block  grid  was  used  where  a  separate  grid  block  extended  from  the  blunt 
trailing  edge  of  the  model  to  the  downstream  boundary.  Far-field  bound¬ 
aries  as  near  as  five  chord  lengths  to  as  far  as  25  chord  lengths  were 
examined  (there  was  no  attempt  to  model  the  wind  tunnel  walls) .  The 
initial  grid  spacing  from  the  surface  ranged  from  0.00005  to  0.001  chord 
lengths.  Riemann  invariant,  non-ref lective  and  free  stream/no-change 
farfield  boundary  conditions  were  all  tried  at  some  point  in  the  study. 
The  net  result  of  these  attempts  was  that,  only  at  Mach  =  0.752  and 

below  could  a  solution  that  matched  the  experimental  data32  be  obtained 
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using  the  basic  equations  and  procedures  described  in  the  previous  sec¬ 
tion.  Figure  3  shows  a  typical  Mach  =  0.752  comparison  of  the  pressure 
coefficients  from  the  MacCormack  explicit  code  used  in  this  study 

(■Shang*  code)  and  the  NASA/Ames  ARC2D  code;  an  implicit,  thin-layer 
procedure.  The  ARC2D  code  is  widely  used  in  government  and  industry  and 
was  used  in  this  study  as  a  "control*  code  for  the  current  study. 

At  M  =  0.807  the  computational  solution  did  not  match  experimental 
data  for  any  of  the  computational  domains  and  boundary  conditions  at¬ 
tempted.  There  were  two  Solutions*  for  the  Mach  =  0.807  condition.  An 
attached  flow  solution,  which  occurred  when  the  eddy  viscosity  was  al¬ 
lowed  to  become  large  (the  attached-flow  equation  for  Fwake  (47) )  ,  and  an 

oscillatory  solution  when  eddy  viscosity  was  kept  small  (the  separated- 
flow/wake  equation  for  F  ,  (48)).  The  oscillatory  solution  was  also 

observed  with  the  ARC2D  code,  which  employs  the  Baldwin-Lomax  turbulence 
model.  Figure  4  shows  how  the  lift  coefficient  varied  with  increasing 
iterations  for  both  ARC2D  and  the  Shang  code.  (The  shape  of  the  oscil¬ 
latory  lift  curve  was  sensitive  to  grids,  boundary  conditions  and  other 
factors.  Comparisons  with  ARC2D  results  are  provided  only  to  show  that 
the  problems  are  not  unique  to  the  Shang  code.) 

Oscillatory  solutions  appeared  to  alternate  between  being  attached 
on  one  surface  and  separated  from  the  other,  then  reversing.  This  ap¬ 
parent  periodic  variation  in  the  lift  coefficient  was  observed  even 
though  the  solutions  were  not  run  time -accurate.  A  similar  periodic 
phenomenon  has  been  observed  experimentally  on  an  18%  thick  circular  arc 
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airfoil,  however,  no  such  oscillations  were  observed  during  the  test  of 
Ref.  32  (according  to  private  conversations  with  the  author).  Lift 
variations  of  the  magnitude  indicated  by  Figure  4  could  have  hardly  gone 
unnoticed. 

Researchers  familiar  with  algebraic  turbulence  models  would  prob¬ 
ably  not  be  surprised  by  the  poor  results  at  Mach  =  0.807.  Indeed, 
every  analysis  of  transonic  turbulent  flow  prediction  consulted  (cf. 
Refs.  16,17,18)  emphasized  the  inadequacies  of  turbulence  models  for 
separated  flow.  This  study  had  to  this  point  only  succeeded  in  rein¬ 
forcing  those  assessments.  Some  success  in  computing  axisymmetric  flow 
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with  local  separation  had  been  reported  by  using  a  "hybrid*  eddy- 

viscosity /Reynolds  shear-stress  turbulence  model.  That  procedure  was 
incorporated  into  the  Shang  code  in  place  of  the  Baldwin-Lomax  model, 
but  the  results  were  not  noticeably  better  than  the  previous  results. 

Of  course  we  must  always  be  skeptical  of  results  obtained  from  trying  to 
duplicate  another  author’s  procedure  with  only  limited  information. 

Analysis  of  the  Algebraic  Turbulence  Model 


Some  insight  into  the  oscillating  solution  problem  was  gained  by 
looking  at  the  results  for  Mach  =  0.75  on  a  340  x  85  grid,  and  with  the 
Baldwin-Lomax  model  as  described  in  Ref.  11.  That  is,  over  the 

body  was  chosen  as  the  minimum  of ; 


wake 


=  (Y  F  ), 
v  max  max7 


or, 


^wake  ^wk  ^max  udif 


0/*, 


max 


The  convergence  proceeded  rather  rapidly;  however,  the  lift  and  drag 
never  reached  a  constant  value  (Figure  5) .  There  was  evidence  of  small 
oscillations,  primarily  in  the  shock  location  on  the  upper  surface 
(Figure  6) .  A  change  was  made  in  the  program  so  that  Pwajte  w&s  defined 

only  as  (Y  F  )  over  the  body  surface,  and  the  oscillations  im- 

v  max  max7  7  ’ 

mediately  disappeared.  Apparently  the  increase  in  dissipation  from  the 
greater  eddy  viscosity  on  the  upper  surface  was  sufficient  to  fix  the 
shock  location  and  permit  convergence. 

The  ineffectiveness  of  changes  in  the  computational  domain,  the 
acknowledged  inadequacies  of  simple  eddy-viscosity  models,  and  the  above 
results  at  M  =  0.75,  suggested  that  the  key  to  converged  solutions  for 
the  Mach  =  0.807  condition  might  lie  in  establishing  the  proper  level  of 
eddy-viscosity  for  the  separated  region  aft  of  the  shocks.  A  more 
detailed  investigation  of  the  key  parameters  in  the  Baldwin-Lomax  model 
was  initiated.  In  the  attached-flow  formulation  (equation  (47)), 

outside  of  the  inner  layer  is  essentially  a  function  of  Yfflax  and  Ivl. 
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P  .  '  Y2  IWIV 

wake  max  Y 

max 

Therefore,  when  the  boundary  layer  is  thin,  Yaay  is  small  and  thus 

is  small.  In  separated  flow  behind  a  shock,  however,  the  boundary  layer 
thickens  rapidly,  Ymav  gets  large,  and  Pwa^e  increases  as  the  square  of 

The  resulting  large  values  for  eddy  viscosity  produce 

large  amounts  of  dissipation  in  the  separated  flow  region  and  suppress 
any  tendency  for  separation.  This  results  in  the  strong  attached-flow 
convergence  for  the  M  =  0.752  and  M  =  0.807  (Figure  8). 

The  separated-f low/wake  formulation  for  Fwake’  (equation  (48) ) , 

outside  of  the  inner  layer  is  essentially; 


Y  (Figure  7) 
max  v  °  1 


wake 


-.k  “di( 


max 


C  .  =  constant. 


Notice  that  boundary -layer  thickness  does  not  appear  in  the  expression. 
Rather,  F^^  varies  as  the  square  of  the  maximum  velocity  at  the  x 

location  of  interest.  Figure  9  shows  t w I y  and  F  ,  for  a  typi- 

max 

cal  oscillatory  result  for  U  =  0.807.  decreases  rapidly  aft  of  the 

shock  (approximately  60  percent  chord) .  Since  F^a^e  is  a  function  of 

the  square  of  changes  in  velocity  would  tend  to  have  a  large  and 

immediate  effect  on  F  ^  ,  and  thus,  eddy-viscosity. 

Consider  the  following  hypothesis:  Some  perturbation  causes  the 
shock  wave  on  one  surface  to  become  slightly  stronger  which  reduces  the 
velocity  downstream  of  the  shock.  The  vorticity,  u,  at  Ymax  may  also 

decrease,  but  since  is  proportional  to  the  square  of  the 

immediate  effect  is  that  eddy-viscosity  is  reduced.  The  resulting 
decrease  in  dissipation  aft  of  the  shock  produces  a  larger  separated 
flow  region,  drives  the  shock  farther  forward,  and  it  becomes  still 
stronger.  This  produces  an  even  lower  and  the  process  continues. 

Eventually  something  (change  in  circulation,  vorticity,  etc.)  stops  the 
forward  movement  and  the  shock  starts  to  move  aft  and  weaken.  Again, 
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the  process  is  self-sustaining  because  increases,  eddy-viscosity 

gets  larger,  and  the  separated  region  decreases.  The  shock  continues  to 
move  aft  until  the  flow  is  fully  attached  on  that  surface,  and  then  the 
cycle  repeats. 

This  type  of  self-sustained,  oscillatory  behavior  was  observed  for 
all  computations  that  had  strong  shocks  with  massive  separation  on  both 
upper  and  lower  surfaces.  However,  it  was  not  obvious  during  the  ex¬ 
periment;  at  least  not  to  the  magnitude  indicated  by  the  computed 
results.  Obviously,  something  was  missing  in  the  computed  results,  and 
all  indications  were  that  the  turbulence  model  was  not  providing  correct 
levels  of  eddy-viscosity.  Something  was  needed  in  the  definition  of 
^wake  defeated  the  self-sustaining  movement  of  the  shock  wave  and 

the  corresponding  change  in  the  separated  flow  region. 

Modifications  to  the  Algebraic  Turbulence  Model 


Examination  of  the  flow  in  the  separated  region  (Figure  10) , 
revealed  a  reversed  flow  region  near  the  body,  and  a  high  velocity 
gradient  region  at  the  outer  edge  of  the  separated  flow.  An  attempt  was 
made  to  treat  these  two  regions  with  separate,  attached-flow  turbulence 
models  by  applying  equation  (47)  to  each  region  independently.  That 
approach  proved  unsuccessful  because  the  eddy-viscosity  became  so  large 
that  the  reversed  flow  region  was  completely  suppressed.  Another  ap¬ 
proach  based  on  the  separated -flow/wake  formulation  (equation  (48))  for 
Fwake  hbe  separated  region,  plus  a  weighted  function  of  the  height  of 

the  reversed  flow  region  to  represent  the  additional  dissipation  needed, 
was  developed.  The  resulting  equation  may  be  written  as; 


F  ,  =  ^wk  ^max  ^dif  +  (D  ,Y.)  ^wk  ^max  ^dif 

wake  - p -  wk  O'  - p - 

max  max 

The  last  term  is  the  additional,  weighted  dissipation,  and  is  propor¬ 
tional  to  the  height  (Y  )  of  the  reversed  flow  region. 


The  original  equation  for  F^^g 


can  be  retained  if  a  new,  variable 


coefficient  is  defined  as; 
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°.ks  *  C„kP.k  \  *  «• 

Where, 

is  the  height  of  the  reverse  flow  region, 
Dwk  is  a  weighting  function  whose  value  is 


to  be  determined. 

Cwks  was  su^s^itute^  f°r  Cwk  in  the  equation  for  Fwake  and  was  computed 
at  each  x  location  for  each  iteration.  The  value  of  C  ^  changes  in 

direct  proportion  to  the  height  of  the  reverse  flow  region,  which 
provides  the  additional  dissipation  necessary  to  counteract  the  in¬ 
fluence  of  changes. 


The  new  parameter,  ,  was  evaluated  to  determine  the  magnitude 

needed,  and  to  assess  its  effect  on  the  solution.  Values  from  20  to  120 
were  used.  Figure  11  shows  the  effect  on  Ct  of  the  different  values  for 


D^.  As  would  be  expected,  the  larger  values  of  provide  large  eddy- 

viscosity  and  produce  results  that  approach  the  solutions  obtained  with 
the  attached  flow  equation  for  The  smaller  values  of  provide 

solutions  that  show  signs  of  the  oscillations  observed  with  the  original 
separated-f low/wake  equation  for  ?wajce •  Figure  11  indicates  that  the 

lift  approaches  a  minimum  value  at  about  Dwk  =  60,  and  oscillates  about 

this  mean  value  at  lower  levels  of  D^*  Unfortunately,  drag  (Figure  12) 

does  not  seem  to  approach  a  fixed  value  as  Dwk  is  reduced.  Rather,  the 

mean  drag  continues  to  drop,  while  the  magnitude  of  the  oscillations 
increase,  at  the  lower  values.  This  is  probably  a  hysteresis  ef¬ 
fect,  aggravated  by  the  extremely  fine  grid  spacing  near  the  surface, 
that  would  vanish  at  a  large  number  of  iterations.  The  pressure  coeffi¬ 
cient  data  is  not  strongly  affected  by  the  oscillations  introduced  by 

low  values  of  D^.  Figure  13  provides  pressure  data  at  four  locations 

of  an  oscillation  when  equals  20.  There  is  very  little  effect  on 

shock  location  or  the  flow  upstream  of  the  shock.  What  effects  there 
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are  appear  to  be  confined  to  the  foot  of  the  shocks  and  the  nearby 
separated  regions. 

New  Results 


As  a  final  test  of  the  modified  separated-f low/wake  eddy-viscosity 

model,  a  new  solution,  starting  with  free-stream  conditions  everywhere, 

♦ 

was  initiated  for  the  Mach  =  0.807,  a  =  0.36,  Reynolds  number  =  6.2  x 


106  case  with  =  SO.  This  value  was  chosen  for  because  it  ap¬ 

pears  to  be  large  enough  to  provide  a  non-oscillatory  solution,  yet 
small  enough  to  be  in  the  constant  range  for  C^.  The  solution  was  es¬ 


sentially  converged  after  20,000  iterations.  Another  20,000  iterations 
were  computed  just  to  make  sure  the  convergence  was  really  stable 
(Figure  14).  Figure  15  shows  pressure  coefficient  results  compared  to 
experimental  pressure  coefficients.  Converged  solutions  were  also  ob- 


6  6 

tained  for  Reynolds  numbers  of  12.0  x  10  and  2.0  x  10  to  evaluate 
Reynolds  number  effects  on  the  value  for  D^.  The  value  for  should 

be  the  smallest  value  that  provides  a  non-oscillatory  (or  acceptably 
small  oscillatory)  solution. 
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III.  CONCLUSIONS 


Governing  equations  for  the  flow  of  an  ideal  gas,  and  a  review  of 
the  process  involved  in  a  numerical  solution  of  the  resulting  Navier- 
Stokes  equations  has  been  provided.  This  includes  a  discussion  of 
finite  difference  algorithms,  boundary  conditions,  computational  mesh 
(grid)  generation  and  the  development  of  an  algebraic  turbulence  model. 

Unfortunately,  numerical  solutions  of  the  governing  equations  for 
flow  fields  with  turbulent  boundary  layers  and  strong  shock-induced 
separation  that  match  experimental  data  has  been  beyond  the  state  of  the 
art.  The  limiting  problem  was  associated  with  the  closure  model  for  the 
Reynolds  averaged  Navier-Stokes  equations.  This  report  describes  a 
modification  to  the  Baldwin-Lomax  algebraic  eddy-viscosity  turbulence 
model  that  appears  to  do  a  very  good  job  of  predicting  shock  location  on 
both  the  upper  and  lower  surfaces.  The  computed  shock  is  not  quite  as 
sharp  as  the  experimental  data,  but  this  could  be  due  to  a  lack  of  grid 
resolution  at  the  shock  location.  The  shape  of  the  pressure  distribu¬ 
tion  in  the  separated  region  is  well  modeled;  however,  the  level  of  the 
pressures  are  somewhat  high  (low  negative  pressure  coefficients).  There 
also  appears  to  be  a  small  discrepancy  in  the  pressures  on  the  upper 
surface  ahead  of  the  shock.  The  boundary  conditions  imposed  by  the  wind 
tunnel  walls  could  have  a  significant  effect.  Also,  grid  resolution  can 
be  contributing  to  the  lack  of  agreement.  The  grid  used  for  this  study 
clusters  points  near  the  surface,  which  would  be  more  appropriate  for  an 
attached,  viscous  flow.  When  the  flow  separates  the  high  gradient 
velocity  region  lifts  off  the  surface  (Figure  16)  and  is  in  an  area  with 
relatively  sparse  grid  points. 

Additional  research  is  needed  to  verify  the  procedure  for  different 
airfoils  and  test  conditions.  Additional  work  is  also  needed  to  inves¬ 
tigate  values  for  two  new  variables  in  the  turbulence  model,  and  Yq . 

At  this  time  the  recommended  approach  for  defining  is  to  pick  the 

smallest  value  that  will  provide  a  solution  with  acceptably  small  oscil¬ 
lations  in  the  shock  location.  The  choice  of  Yo  as  the  height  of  the 
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reverse  flow  region  is  also  rather  arbitrary.  Any  length  scale  that  is 
proportional  to  the  magnitude  of  separation  could  be  used.  Perhaps  with 
some  additional  study  a  combination  of  Yo  and  D  ^  that  would  work  for  a 

wide  variety  of  airfoils  could  be  determined.  The  applicability  of  this 
procedure  to  three-dimensional  flow  also  needs  to  be  investigated. 
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Table  1 


Supersonic 
Flow  at 
Boundary 


Subsonic 
Flow  at 
Boundary 


*  Courtesy 


Riemann  invariant  Far  Field  Boundary  Conditions 


* 


In  Flow  Boundary  Out  Flow  Boundary 

-  |  - 

All  variables  set  to  free  I  All  variables  extrapolated 

stream  values.  I  from  the  interior. 


V.  . .  ,  and  entropy  I  V,  ...  and  entropy 

tangential  I  tangential 

set  to  freestream  values.  I  extrapolated  from  interior. 

^normal  ant*  sPee^  SOUI,d  computed  from  fixed  and  extrapolated 
Riemann  invariants. 

Density  and  pressure  calculated  from  entropy  and  speed  of 
sound . 


of  Pradeep  Raj,  Lockheed-California  Company 
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Table  2 


Values  for  Parameters  in  Grid  Code 


Variable 

Computer 

Symbol 

Value 

Comment 

i 

max 

JMAX 

340 

Number  of  points  around  the 
body  and  wake  outline. 

^max 

KMAX 

85 

Number  of  grid  levels  away 
from  the  body  to  the  far  field. 

hr) 

SETA 

0.0001 

Initial  step  size  away  from 
the  body. 

i) 

'max 

SETMX 

15 

Far  field  location  in  chord 
lengths. 

ESCAL 

0.00127 

Factor  controlling  the  rate  at 
which  the  grid  attempts  to  pro¬ 
vide  equal  cell  areas. 

SMUIM 

0.005 

Explicit  smoothing  parameter. 

ai 

SMU 

0.005 

Implicit  smoothing  parameter. 

a 

ALPHA 

1  to  2.5 

Factor  that  changes  the  finite 
difference  approximation. 

Varies  as  a  function  of  grid 
level . 

I 
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C-Type  Grid 


0-Type  Grid 


Figure  1.  0-Type  and  C-Type  Grids  on  NLR-7301  Airfoil 


MACH  CONTOURS 


Figure  2.  NLR  7301  Airfoil  with  Strong  Shocks  and 

Separated  Flow,  Mach=0.807,  a=0.37,  R  =12X10 


SURFACE  PRESSURE  COCrriCIENTS 


Shang  Code  Solution  ARC2D  Solution 


ITERATIONS  »  ITERATIONS 


Figure  4  Oscillating  Solutions  at  Mach  =  0.807 


ITERATIONS 


Figure  5  Oscillations  at  Mach=0.75 


Mach=0.75  Range  of  Pressure  Distributions 


FWAKE  COMPARISON 


YMAX  LOCATION 


Figure  7  and  Ymax  for  Mach=0.807,  a=0.37,  Rn=6.2X106 
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V0RT1CITY 


KE  COMPARISON 


•  *1 


U  M  14 

X/C 


for  Mach=0.807  Oscillatory  Flow 


wake  velocity  profile 


WAKE  VELOCITY  PROFILE 


>•■9  tov  toe  BJ  til  1*9  «&7  MB  HI  M.3  309 

ITERATIONS  mlCf 


Figure  11  Effects  of  on  Lift,  Mach=0.807,  a~ 0.37,  Rn=6.2X10 


DRAG  COEFFICIENT  VS  ITERATIONS 


Figure  12  Effects  of  on  Drag,  Mach=0.807,  a=0.37,  Rn=6.2X106 
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LIFT  COEFFICIENT 


PRESSURE  OSCILLATIONS 


Experimental  Data 
A  Upper  Surface 
A  Lower  Surface 


Computed  Results 


LIST  OF  SYMBOLS 


A 

B 

C 

C 


cp 


"kleb 


Jwk 


wk 


kleb 


max 


wake 


I 

f 

G 

i 

j 

n 

K 

£ 

k 

k. 


Surface  area  of  an  arbitrary  volume 

Sublayer  constant  for  turbulence  model 
Matrix  defined  in  grid  development 
Matrix  defined  in  grid  development 
Constant  in  turbulence  modeling 

Klebanoff  constant  in  turbulence  modeling 

Lift  coefficient 

Constant  pressure  specific  heat 

Constant  volume  specific  heat 

Constant  in  turbulence  modeling 

Speed  of  sound 

van  Driest  damping  function 

Variable  in  turbulence  modeling 

Flux  vector  in  governing  equations 
Klebanoff  intermittency  function 

Function  in  turbulence  modeling 

Function  in  turbulence  modeling 

Right  hand  side  vector  in  grid  development 

Arbitrary  variable 

Flux  vector  in  governing  equations 

Index  representing  current  x  or  (  location 

Index  representing  current  y  or  1]  location 

Index  representing  current  time  step 

Clauser  constant  in  turbulence  modeling 

Thermal  conductivity  coefficient 

von  Karman  constant 

Eddy  conductivity  coefficient 
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t 

p 

Prt 

$ 

•¥ 

q 

i 

R 

R 

n 

T 

t 

0 

U 

u 

V 

v 

X 

y 

Greek: 


Dyadic  representing  stresses 
Hydrostatic  pressure 
Turbulent  Prandtl  number 

Arbitrary  tensor 

Heat  flux  vector 

Vector  in  grid  development 
Gas  constant 
Reynolds  number 

Temperature 
Time  dimension 

Velocity  vector 

Vector  of  flow  variables  in  governing  equations 

Velocity  in  the  x  direction 

Arbitrary  control  volume 

Velocity  in  the  y  direction 

Cartesian  coordinate 

Cartesian  coordinate 


a  Constant  in  MacCormack  pressure  damping 

a  Angle  of  attack 

6  Finite  difference  operator 

e  Eddy-viscosity  coefficient 

rj  Transformed  coordinate  for  y 

7  Ratio  of  specific  heats 

t  Length  scale  for  turbulence  modeling 

X  Coefficient  of  bulk  viscosity 

H  Coefficient  of  shear  viscosity 

/)  Density  (waoo  pnr  unit  volume) 

o  Smoothing  parameter  in  grid  development 

r  Stress  on  a  fluid  element 


64 


v  Vorticity 

4  Transformed  coordinate  for  x 

A  Incremental  values 

A  Dissipation  operator  in  grid  development 

V  Divergence  operator 

V  Dissipation  operator  in  grid  development 
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